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We consider a microfluidic mixer based on hydrodynamic focusing, which is used to initiate 
the folding process of individual proteins. The folding process is initiated by quickly diluting 
a local denaturant concentration, and we define mixing time as the time advecting proteins 
experience a specified to achieve a local drop in denaturant concentration. In previous work, 
we presented a minimization of mixing time which considered optimal geometry and flow 
conditions, and achieved a design with a predicted mixing time of 0.10 ps. The aim of the 
current paper is twofold. First, we explore the sensitivity of mixing time to key geometric 
and flow parameters. In particular, we study the angle between inlets, the shape of the 
channel intersections, channel widths, mixer depth, mixer symmetry, inlet velocities, working 
fluid physical properties, and denaturant concentration thresholds. Second, we analyze the 
uniformity of mixing times as a function of inlet flow streamlines. We find the shape of 
the intersection, channel width, inlet velocity ratio, and asymmetries have strong effects on 
mixing time; while inlet angles, mixer depth, fluid properties, and concentration thresholds 
have weaker effects. Also, the uniformity of the mixing time is preserved for most of the 
inlet flow and distances of down to within about 0.4 pm of the mixer wall. We offer these 
analyses of sensitivities to imperfections in mixer geometry and flow conditions as a guide 
to experimental efforts which aim to fabricate and use these types of mixers. Our study also 
highlights key issues and provides a guide to the optimization and practical design of other 
microfluidic devices dependent on both geometry and flow conditions. 

Keywords: Microfluidic mixers; Device design; Numerical Modeling; Sensitivity analysis; Mixing 
uniformity. 
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I. INTRODUCTION 


Protein folding studies [10] pose a great challenge to microfluidic mixers. Proteins are composed 
of chains of amino acids which take on complex three-dimensional (3D) structures to achieve a wide 

riQ 

range of biological functions [1|, yj. The range of applications of protein folding in research and 
industry is wide and includes drug discovery, DNA sequencing, and molecular analysis or food 
engineering [s-^. One of the most versatile methods of initiating the process of protein folding is 
using changes in chemical potential (e.g., changing the concentration of a chemical species) 0-l8(. 
In this work, we consider a class of microfluidic mixers based on diffusion from (or to) a hy- 

u 

drodynamically focused stream. This type of mixer was initially proposed by Brody et al. |0. A 
geometrical representation of such a mixer is shown in Figure [2] (here our optimized mixer of Ref. 
171 ). The basic features of the design are as follows: It is composed of three inlet channels and 
a common outlet channel, and the geometry has a symmetry with the center channel. Typically, 
a mixture of unfolded proteins and a chemical denaturant solution is injected through the center 
channel and exposed to background buffers (no denaturant) streams through the two side chan¬ 
nels. The design goal is to rapidly decrease the denaturant concentration in order to rapidly initiate 

plication of Brody et ah, there have been sig- 


protein folding in the outlet channel [10|| . Since the pu 
nificant advances on the design of these mixers 


Ill4l3l| including reduction in consumption rate of 


reactants, methods of detection, manufacturing and, perhaps most importantly, drastic reductions 
of the mixing time (i.e. the time required to reach a sufficiently low denaturant concentration). 

We recently studied the optimization of the shape and flow conditions of a particular hydro- 
dynamic focused microfluidic mixer [l^. The objective was to improve the mixing time of the 
best mixer designs found in literature, which exhibited mixing times of approximately 1.0 ^s. To 
this end, we introduced a mathematical model which computes the mixing time for a given mixer 
geometry and injection velocities. Then, we defined the corresponding optimization problem and 
solved it by considering a hybrid global optimization method 1J-[10|. This approach was carried 


out and presented using both 2D and 3D models. To save on computational time, much of the 
optimization process was conducted using the 2D model. However, our earlier work also pointed 
out that certain important effects (including the impact of upper and lower mixer walls and inertial 
effects on the velocity field) can be appreciated only with the 3D model. We therefore performed 
3D model studies to analyze such effects. The optimized mixer generated by our approach achieved 
a mixing time of about 0.10 fis. The shape of this optimized microfluidic mixer, its concentration 
distribution and the concentration evolution of a particle in its central streamline are summarized 
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FIG. 1. Optimized mixer of Reference Q and associated mixing performance: (a) top view representation of 
the half of the mixer shape (symmetry with respect to x=0 /im ) of the mixer with a superposed grey scale 
plot of the denaturant normalized concentration distribution c at width z=0 ^m and (b) the time evolution 
of the denaturant concentration of a particle flowing along the symmetry streamline starting from x=0 /rm, 
y=5 /im and z=0 /im. The parts of the mixer shape corresponding to the protuberances Protup and Protjo, 
introduced in Section PV A 21 are highlighted in sub-figure (a). 


in Figure [TJ The optimized side and center channel injection velocities were Ug =5.2 m s“^ and 
Uc =0.2 m s“^, respectively. The optimization problem studied in this previous work was identihed 
as highly nonlinear 2^. Further, the process has many parameters which are difficult to know with 


great precision in experiments. Therefore, it is important to understand and quantify the stability 
of the device performance with respect to these parameters. This enables identification of key 
parameters and so guide experimental efforts. To this aim, we have analyzed the optimized mixer 
to study and quantify its robustness to parameters perturbations. 

We previously presented a very simple sensitivity analysis [l^. That preliminary sensitivity 
study consisted of random perturbations of all the parameters by taking uniform variations within 
a range of [—/3%, +/3%] of their value. Results showed that the mixing time variations were of 
the same order as the normalized perturbations considered, suggesting the optimized solution was 
fairly stable. See Ivorra et al [l^ for further information. 
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Here, we significantly increase the scope of our sensitivity analysis. We quantify the impact of 
mixing time on the key design parameters of the mixer. The objective of our study is to provide 
recommendations and guidelines for the fabrication of the device introduced here. More precisely, 
we consider and study (i) Geometrical parameters defining the mixer shape: the angle defined by 
the channel intersection, the shape of the channel intersection, the width of the inlet and outlet 
channels, the mixer depth and possible irregularities in the symmetry of the shape;(ii) Central and 
side injections velocities; and (iii) Physical coefficients associated with the working fluid and the 
concentration thresholds of the mixing time definition. 

In addition to those sensitivity analysis experiments, we also analyze the uniformity of the 
mixing time as a function of the inlet streamline location in the inlet channel. This mixing time 
uniformity analysis quantifies the robustness of the mixing time through the whole inlet flow, and 
helps place a statistical confidence on observed mixing times. In particular, it helps quantify the 
so-called wall effect (due to the no-slip condition at the mixer walls, resulting in low velocity values 
near the walls) on mixer performance. 

The last two decades have seen a large number of microfluidic device designs and their use in 
a wide range of applications. Most, if not all, of these devices have performance specifications 
which are dependent on their geometry and flow control conditions (e.g., flow rates, pressure, inlet 
concentrations). Despite this, the systematic study of how performance depends on intentional or 
untintential design parameters is rarely if ever demonstrated. For this reason, we also offer the 
current work as a case study describing the significant challenge and complexity of determining 
design robustness for microfluidics. 

This article is organized as follows: Section |TI] introduces the 3D model used to estimate the 
mixing times. Section IIIII describes the mixing time uniformity analysis and the results. Section 
IIVI presents the numerical experiments carried out to perform the extended sensitivity analysis and 
deduce major conclusions and design guidelines. 


II. MICROFLUIDIC MIXER MODELING 

We consider the microfluidic mixer described in Section HI The geometry has two symmetry 
planes which we use to reduce the simulation domain to a quarter of the mixer. This reduced 
domain is denoted by D, as depicted in Figure [2j The mixer shape is composed of interpolated 
surfaces, and the inlet velocities are described by a set of parameters denoted by 4>, detailed in Ref. 
17. 
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Side Inlet 



FIG. 2. Typical three-dimensional representation of the microfluidic mixer geometry. The mixer design 
hydrodynamically focuses a center inlet stream using two side inlets. In dark gray we represent the domain 
used for numerical simulations. The geometry’s two symmetry planes are also highlighted. 


We consider guanidine hydrochloride (GdCl) as the denaturant 


liquid flow is incompressible 


Q,Ei. 


We assume the mixer 


12l |. Thus, the flow velocity and the denaturant concentration distri¬ 


bution are approximated by using the steady configurations of the incompressible Navier-Stokes 
ed with the convective diffusion equation. More precisely, we consider the following 

3 : 


equations coup 


system [ll|, [l^. 


—V • (r/(Vu -h (Vu) ) — pi) -h /9(u • V)u = 0 in Id, 


V • u = 0 


in n, 
in n. 


( 1 ) 


V • {—D'Ve) -h u • Vc = 0 
where c is the denaturant normalized concentration distribution, u is the flow velocity vector (m 
s“^), p is the pressure held (Pa), D = 2x 10“® is the diffusion coefficient of the denaturant solution 
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in the background buffer (m^ s“^), ry = 9.8 x 10““^ is the denaturant solution dynamic viscosity (kg 
s“^) and p = 1010 is the denaturant solution density (kg ). 

System ([T]) is completed by the following boundary conditions: 

For the flow velocity u: 


u = 0 on 

u = —ttsparain on F*, 

< u = —ttcpara 2 n on Fc, 

p = 0 and (ry(Vu + (Vu)''"))n = 0 on Fg, 

n • u = 0 and t • (ry(Vu + (Vu)"*") — pl)n = 0 on Fq, 


( 2 ) 


where Fc, F^, Fg, F^ and Fq denote the boundaries representing the central inlet, the side inlet, 
the outlet, the mixer walls and the symmetry plane, respectively; Ug and Uc are the maximum side 
and center channel injection velocities (m s“^), respectively; parai and para 2 are the laminar flow 
profiles, which are equal to 0 in the inlet border and to 1 in the inlet center, of the side and central 


inlets, respectively 


19|; and (t, n) is the local orthonormal reference frame along the boundary. 


For the concentration c: 


n • {—D'Ve + cu) = —cqu on Fg, 
c = 0 on Fg, 

< 

n • {—DVe) = 0 on Fg, 

n • {—DVe + cu) =0 on F^ U Fq, 


( 3 ) 


where cq = 1 is the initial denaturant normalized concentration in the center inlet. 

In this work, the mixing time of a particular mixer (j), denoted by J{4>) is defined as the time 
required to change the denaturant normalized concentration of a typical Lagrangian stream fluid 
particle situated in the symmetry streamline at depth z = 0 pm from a% to uj%. It is computed 
by: 


m 



dy 

u'^(y)’ 


( 4 ) 


where u'^ and denote the solution of System dU-®, when considering the mixer defined by 
(j)‘, and Co and ct denote the y-coordinate of points situated along the streamline defined by the 
intersection of the two symmetry planes z = 0 pm and x = 0 pm, i.e. the y-axis, where the 
denaturant normalized concentration is a and oj, respectively. By default, we assume a = 90% 
and uj = 30%. 
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The numerical model used to approximate the solutions of System ([I])-© and to compute dH) 
was implemented by coupling Matlab scripts with COMSOL Multiphysics 3.5a models. 

III. UNIFORMITY OF THE MIXING TIME 

We first analyze the non-uniformity of mixing times across the focused stream for our optimized 
mixer (j)o. Indeed, as suggested in Refs. 13 and ll^. the mixing time can be measured not only in 
the symmetry streamline, situated on the (x,z)=(0,0) segment, but also in other streamlines. We 
are interested in the uniformity of mixing times, as protein states in these mixers are quantified 
experimentally within a finite probe volume which integrates signal throughout a volume in space 
within the mixing region. This measurement volume is fed in principle by all streamlines of the 
center inlet channel. 

We consider 100 streamlines, denoted by starting from a finite set of points, which 

are denoted by in Tc. Here, = 1,...,10 and j = 1,...,10} where P(ij) = 

(j^0.9/rm, |0.75/im). In the previous definition, the maximum coordinate in the x-axis (i.e., 0.9 
/im) has been selected in order to avoid particles too close to the wall T^, and the maximum 
coordinate in the z-axis (i.e., 0.75 /rm) has been chosen as a characteristic 1.5 /rm depth of field 
for confocal microscope imaging (i.e., extent of the measurement volume) Ut]. Those streamlines 
are numerically approximated by considering an explicit Euler scheme and the velocity vector u 
obtained by solving System ®-o B- 

For each streamline sltj, we compute the associated mixing times, denoted by in a manner 
similar to Equation dH) . More precisely, . j is defined as the time required by a protein within a 
Lagrangian fluid particle to travel from Cqq’^ to c^q’^ , where Cqq’^ and c^q’^ denote the points within 
slij with a concentration of 90% and 30%, respectively. Next, we study the spatial distribution 
according to the streamline starting point in the maximum value, the mean value and the 
standard deviation of {tsk j)i^j=i- Furthermore, we also compute the weighted mixing time value 


of slij, denoted by tgi^ and defined as 


t q]■ . - 


2 ^i,j = l ^i,j 


(5) 


where Wjj denotes the velocity of a particle in the streamline slij at its initial position This 

choice of weight coefficients reflects the fact that the probe volume used to measure experimentally 
the mixing time receive particles more frequently from streamlines with the highest velocities. 


The maximum and standard deviation values of those weighted mixing times {tsii j)ij=i are also 
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studied. 


Furthermore, due to the fact that the depth of the mixer is 10 times larger than the minimum 
width of the center channel, the mixing time variations in the z-axis direction are negligible in 
comparison to the variations in the x-axis [ll]. Thus, we perform a more extensive uniformity 
analysis along the x-axis, by considering 100 streamlines, denoted by (s/i^ 2 =o)i£i, ™ plane 
z = 0 starting from the set of points P(ij) = (p^0.9/im, O/rm) in Fc- The methodology is the same 
as that introduced previously. In this case, we also compute the evolution of both the mean value 


and standard deviation of {tsU and (t.s/.-, with k = 1,...,100. These results will be 

compared with the ones presented in Ref. 111. 


Our study of mixing time uniformity yielded that the mean mixing time value obtained by 
considering {tgi^ was 0.34 //s with a standard deviation of 0.17 /xs. As expected, the maximum 

mixing time value was reached at the streamline s/io,io with a value of 1.43 /xs. 


The mixing times and weighted mixing times of the considered streamlines 
are presented in Figure [3]-( a). As shown, within 0.4 /xm of the centerline, the mixing times vary 
between 0.1 /US and 0.5 /xs, and this region accounts for 60% of the detection events (i.e., considering 
the sum of the weight coefficients (X]i=i j)/j=i contrast, the near-wall region 

of [0.7, 0.9] /xm of the centerline have mixing times between 1 and 1.43 /xs, but these streamlines 
contribute to only 10% of detection events (i.e., considering 


Next, the mixing times tgi. and weighted mixing times tgii across the streamlines 
(sZj 2 =o)j^£i plotted versus spanwise streamline position in Figure [3l-(b). For these 100 stream¬ 
lines, the mean mixing time computed by considering with a standard 

deviation of 0.16 /xs. Again, we can observe that particles near the walls exhibit higher mixing 
times (>1 /xs). However, these near-wall-slow-moving particles contribute only infrequently to 
probe volume detection events. 


We note similar phenomena were reported in Ref. [llj. However, the mixer presented in that 
work exhibited a mean mixing time, considering streamlines in the plane z=0, of 3.1 /xs with 
a standard deviation of 1.5 /xs. The maximum mixing time value was 10 /xs, obtained for the 
streamline closer to the wall F^. The optimized mixer design presented here therefore offers better 
mixing time uniformity leading to more consistent measurements and less scatter in measurement 
ensembles. 
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FIG. 3. Results obtained during the analysis of the mixing time non-uniformity for the optimized mixer: 
a) Mixing and weighted mixing times obtained according to the position (x, z) in Fc of the initial particle 
for the considered streamlines, b) Mixing and weighted mixing times obtained as a function of the position 
X in Fc of the initial particle and for the streamlines considered in the plane z = 0. The weighted mixing 
time reflects the frequency of events (measurements of proteins along said streamline) as determined by 
the stream-line averaged velocity. The lower velocities near the wall yield longer mixing times but are less 
frequent. 
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IV. SENSITIVITY ANALYSIS OF THE MODEL PARAMETERS 

We here present a study of the influence of key parameters of the model described in Section 
m on mixer mixing time. We vary parameters individually, fixing the values of others to the 
corresponding value of the optimized mixer (po- We note that, in our previous work, we explored 
the impact of simultaneous perturbations on the whole set of parameters on mixer performance [l^ . 
We here perform the more complete influence of individual perturbations on the mixing time. 
We believe such individual parameter perturbation analyses are also more useful to designers in 
identifying key parameters and methods for fabrication. We consider the following percent variation 
function: 

— Jm —■ 

where cpp represents the perturbed mixer. 

The parameters analyzed can be classified in three categories: (i) geometrical parameters defin¬ 
ing the mixer shape; (ii) central and side injections velocities; and (iii) physical coefficients associ¬ 
ated with the denaturant solution and the concentration threshold in the mixing time definition. 


A. Geometrical parameters 

In the following computational experiments, we analyze the variation on the mixing time due 
to changes in: (i) the angle defined at the channel intersection; (ii) the shape of the channel inter¬ 
section; (iii) the width of the inlet and outlet channels; (iv) the mixer depth; and (v) perturbation 
in the symmetry of the mixer shape. 


1. Inlet intersection angle 

First, we study the angle between the x-axis and the mixer side channel, denoted by 9. The 
optimized value 9 = vr/S is varied from 0 up to 27r/5 by considering 50 equally spaced interme¬ 
diate values (i.e, we perform 50 evaluations of our model). A geometrical representation of those 
variations is showed in Figure [H 

Perturbations on 9 have generated a mean variation in the mixing time of 3%. Figure [5] gives 
a graphical representation of the obtained results. As we can observe on this plot, the maximum 
variation was around 15% and was obtained for 9 = ‘I'kI'o. Furthermore, the variation was less 
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x-axis (jum) 


FIG. 4. Three mixer shapes for inlets intersection angles of 6* = 0 (light grey), 9 = 7r/5 (white) and 6 = 27r/5 
(dark grey). Only the area where the shape changes is shown. 

than 4% for angles lower than tt/S, and grew up exponentially after that value. This suggests that 
the angle is not a sensible parameter for the mixer performance. 

2. Shape of the ehannel intersection 

We now study the impact of the shape of the area where the three inlets and the outlet intersect. 
The shapes allowed by our model are built by considering Beziers curves and describe a ’bubble’ 
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FIG. 5. Percent variation of mixing time as a function of deviation from the optimal angle 9 value (denoted 
by X=7r/5, c.f. Figure S]). Mixing time is relatively insensitive to small errors in angle of the side channel. 
The distribution shows the strongly non-linear dependence of mixing time on geometry. 

(also called protuberance) invading the central and side inlets from the upper corner (according 
to y-axis) and a protuberance invading the outlet and side inlets from the lower corner. These 
protuberances are dehned according to a restriction (due to a convenient lithographic and plasma 
etching limitation) of a minimum channel width of 1/rm. For the sake of simplicity, those bubbles 
are only described by two scalar numbers Protup and Protio in [0,1], where 0 corresponds to 
the minimum bubble shape and unity is the maximum bubble shape of the upper and lower 
corner, respectively, as allowed by the model. The optimal shape corresponds to Protup=0.8 and 
Protio=0.7. The parts of the mixer shape corresponding to Protup and Protio are presented in 
Figure [TJ A geometrical representation of the minimum, maximum and optimal shapes of the 
protuberances is given in Figure [6j 

This experiment consisted of computing the mixing time of the mixer generated by considering 
all the possible combination of values of Protup and Protio ia [0; 1] with a grid step size of 0.1. This 
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x-axis (iim) x-axis (um) 


FIG. 6. Mixer shapes highlighting the range of mixer shapes we explored. In both cases, the optimal shape 
is shown as a dark line. Shown are detailed views of the shape of the channels intersection near the a) upper 
(denoted by Protup) and b) lower (denoted by Protio) corners of the intersection region. The dark 
gray zones correspond to the domain, between the shape of the maximum and the minimum protuberance 
allowed by the model considered here (according to a minimum channel width of l^m). 


required 121 evaluations of our model. The variation of the mixing time according to analyzed 
values of Protup Protio is presented in Figure [7] and values are reported in Table [H As shown 
by both the figure and the table, for values of Protup and Protio lower than 0.5, the mixing time 
dramatically increased from 50% up to 250%. This indicates that a minimum protuberance in 
both upper and lower corners should be considered in order to obtain an efficient mixing time. 
Furthermore, when Protup and Protio were greater than 0.5, the variation in mixing time was 
moderated and was lowered by 22%, which can be considered as a reasonable value. In addition 
to those first results, we see that the impact on the mixing of Protup was greater than Protio. 
For instance, by decreasing the parameter Protio from 1 to 0 and fixing the value of Protup=l, 
we have generated mixing time variations up to 50%, whereas by decreasing the parameter Protup 
and fixing Protio=l we have obtained a maximum 20% variation of mixing time. This result is 
consistent with the fact that the length of the lower corner is much larger than that of the upper 
(see Figure [6]), thus, its influence on the mixing time is expected to be greater. 

From the previous results, we conclude that the mixing time is sensitive to the shape of these 


protuberances. 
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FIG. 7. Percent variation of mixing time for the optimal shape (represented by X) for the protuberance 
magnitudes considered and described in Section IIV A 21 Protuberance parameter values Protup and Protio 
each vary from 0 to unity with a grid step size of 0.1. The details of the protuberance shape of the side 
channels is an important feature. 


3. Channel width 


We are here interested in estimating the impact of the inlets and outlet widths on the mixing 
time (i.e., the minimum width of these channels where the flow they carry first interact with the 
neighbouring streams). This study is interesting as the mixer design and general shape can be 
scaled geometrically and inserted into different devices. We note that the channel widths were 
fixed during the optimization process in Ref. [ITI and were set to values suited for the mixer 


implementation and validation studies, as the one carried out in Ref. [ll. 


We considered a width denoted by Wc € [l//m,4^m] for the central inlet, a width denoted by Wg € 
[l/xm,4/xm] for the side inlets and a width denoted by Wq E [2/xm,18/im] for the outlet. The original 
optimized shape exhibited Wc = 2/im, Wg = 3/im and Wq = lO^m. All possible configurations 
of channel widths were tested by considering a mesh of step size of 1 /rm for each width, which 
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TABLE 1. Percentage variation in mixing time of the optimal design value and considering the protuberances 
described in Section flV A 21 Here protuberances parameters Protup and Protio, each vary from 0 to unity 
with a grid step size of 0.1. The optimal shape (-) is obtained with Protup=0.8 and Protio=0.7. 


^\Protup 

Protio 

O 

O 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 0.7 0.8 0.9 

1.0 

0.0 

251 220 192 

167 144 126 
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91 

79 
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89 
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218 
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167 145 
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90 

76 

65 

56 

75 

0.2 

186 
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123 
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89 

75 

62 

51 

42 

61 

0.3 

155 
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86 

73 

60 
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39 
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0.4 
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3 
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11 
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represents a total of 272 evaluations of our model. Representations of the mixer shape with all 
channel widths set to their maximum or minimum values, are depicted by Figure [HI 

In order to check the importance of each channel width on the mixing time regarding all possible 
configurations of other width, we considered percent variations denoted hy WE and the mean 
evolution of the mixing time MET according to each width. Both processes are explained below. 
We illustrate the process of computing WE and MET in the case of Wc- This approach can be 
extended to Wg and Wq- 

The value WE^^{j, k) represents a measure of the variation of the mixer mixing time according 
to changes in Wc when other widths are hxed to Wg = j and Wq = k, and is given by 


WE^Sj,k) = 100- 


1 \ J{(j)i,j,k) -meaniJ{4>ij^k)\ 


i=l 




(7) 


where 4 >i,j,k denotes the mixer obtained by considering Wc = i,Ws = j, Wq = k and the other 
parameters set to the optimal values and meauj denotes the mean value of the mixing time 

obtained by varying only i. We compute WEw^{j,k) for j = 1,..,4, k = 2,..., 18 and report its 
mean, minimum and maximum values according to j and k. Those results are reported in Table 

m 
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FIG. 8. Shape of the mixer for lengths of the channels set to their maximum (continuous line for x >0 m) 
values and minimum values (continuous line for a; <0 m) . The dash-dot line corresponds to the axis a; = 0 
(and the y-z symmetry plane). The optimal mixer shape is represented by a dashed line. We note the choice 
of parameterization of the mixer shape (including side inlet channel width) determines the position of the 
region corresponding to the channel intersection. Our geometry variations therefore considered a wide range 
of shapes and relative channel lengths. 


The value MTEw^{i) corresponds to the mean values of the mixing times J{4>ij^k) obtained 
when considering j = 1, ..,4 and k = 2 ,.., 18. The evolution of MTE^^, MTE^^ and MTE^^ are 
depicted in Figure [H 
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FIG. 9. Dependence of mean values of mixing times as defined in Section flV A 31 as a function of inlet and 
outlet channel widths. Shown are plots a) b) MTE^,^ and c) MTE^^ which correspond to the 

mean mixing times obtained when fixing the value of Wc, Ws and Wo, respectively, and let the other widths 
vary. Width variations have a moderate effect on mixing times. 

As we can observe in Table m the most sensitive widths are the side inlets and central inlet 
with a mean mixing time variation of about 65%. This result is expected, since those inlets 
carry the denaturant solution and buffer flows, and thus affect the amount of injected products. 
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TABLE 11. Maximum, minimum and mean values of the mixing time percent variation named WE, defined 
in Section llV A 31 for the widths Wc, Wg and Wo- 



Significant changes to these inlet geometries should be accompanied by changes in inlet velocities 
and performing a new optimization process as in Ref. Il7|. For example, we hypothesize that 
variations which aim to preserve flow rate ratios should be explored first. On the other hand, 
outlet widths in the interval [2,13] //m (from the optimal value of 6 //m) will affect mixing time 
variation by only 7%. Hence, we conclude that such errors on width have only a slight to moderate 
effect on mixing times. Furthermore, regarding Figure[9l we see that the mean mixing time is lower 
when considering values of Wc and Ws in the interval [2//m,4^m] and Wq G [l^m,12/im]. Moreover, 
we remark that configurations with smaller inlets and bigger outlet are the worst from an efficiency 
point of view. 


4- Mixer depth 

Next, we analyzed the effects of the mixer depth (in Z-direction). Imperfections in micro¬ 
fabrication of these mixers can result in depth variations of approximately ztl/rm [l^ . We thus 
computed the mixing time for mixers generated by considering the set of parameter (po and depths 
of 8,9,11 and 12 /rm. The resulting mixing times (and their associated percent variation regarding 
the mixing time of the original mixer with a depth of 10 jim) were 0.14 pLS (34%), 0.12 jjis (13%), 
0.10 /is (6%), and 0.09 /is (13%), respectively. 

As shown by these results, perturbations of ±1 /rm generate reasonable percent variations in 
the weighted mean mixing time between 6% and 13%. As described previously, this indicates that 
errors in the mixer depth due to manufacturing processes do not strongly affect mixing performance 
for these relatively deep (~10 /um) mixers. Note that the highest channel depth yields the lowest 
mixing time (0.09 fis for a depth equal to 12 /rm versus 0.14 //s for the 8 /rm depth). This result 
is expected, as the so-called wall effect (i.e., where the no-slip condition at the top wall results 
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in low velocity values near the wall and near the corner where the X-Y plane meets the Y-Z 
plane) reduces the mixing performance near the mixer walls. Again, we see that the optimal 
mixer design (minimum mixing time) is influenced strongly by changes in manufacturing process 
(namely in achieving high aspect ratio features with deep reactive ion etching). Mixer designs 
with relatively high channel-depth-to-feature width ratios yield optimal results. In our study, the 
minimum channel width (near y = -1.5 /rm) was 1.1 /rm. 


5. Shape Symmetry 

The last geometrical aspect analyzed during this work is the impact of perturbations in the 
symmetry of the mixer according to the plane x = 0 (including nonsymmetric injections velocities) 
on the mixer characteristics. 

To this end, we considered the right half (versus quarter) of the geometry. We then randomly 
generated 100 nonsymmetric mixers by considering perturbations of the parameters from 0.5% 
up to 50% of the left side (respecting to x = 0) of the mixer shape and by keeping the right 
side of the mixer shape to its optimal value. These mixers were then classified according to the 
deviation observed between the streamlines starting from (x = 0, y = 0 ,2 = 0)^m of the symmetric 
and nonsymmetric mixers at the time when the non perturbed symmetric streamline reach Tg. 
According to this classification, we then computed the mean mixing time for each category and 
compared it to the optimized mixer mixing time by considering the percent variation formula ([6]). 
Deviations in the intervals [0,0.3] ym, [0.3,0.6] ym, [0.6,0.9] /rm, [0.9,1.2] ym, [1.2,1.5] ym and 
greater than 1.5 ym generated mean mixing time percent variation of 14%, 64%, 114%, 237%, 
328% and 542%, respectively. 

As we can observe from those data, for deviations below 0.3 /rm, which correspond to parameter 
perturbations lower than 10% in the symmetry of shape and injection velocities, the order of the 
mixing time was conserved with a mixing time variation of 14%. For greater deviations, the mixing 
time was dramatically increased from 64% up to 500%. Thus, we recommend normalized symmetry 
errors of less than 10% be achieved to ensure a mixing time close to the optimal value. 


B. Flow injection velocities 


We studied the influence of injection velocities on mixing time. The optimized injection velocities 
obtained in Ref llTI were Ug =5.2 m s“^ and Uc =0.2 m s“^ (equivalent to a ratio rtc/tts=0.0389). 




20 


For this, we considered the optimized mixer and varied its side injection velocity Ug = from 0.5 
m/s to 9.5 m/s, with a step size of 0.5 m/s. Then, we chose Ug in order to achieve Uc/ug ratios in 
the set {25, 50, 75,100, 250}% (considered as typical values) of the optimal ratio. This part of our 
study required a total of 95 evaluations of our model. 

The results are summarized in Table m and Figure [TOl We can see that for Ug € [4,9] m/s and 
Uc G [0.12,0.88] m/s (i.e., a ratio Ucjug of [0.0292,0.0973]), the mixing time has exhibited variations 
lower than 10%. This suggests our mixer should be robust to small perturbations in the injection 
velocities. In particular, the velocity of the central inlet flow should be in the interval [0.1,0.8]m/s 
to obtain a reasonable mixing time. Moreover, from those results we can deduce that if the ratio 
and/or Ug are too small, the mixing time is drastically increased (more than 1000%). In fact, the 
mixer performance becomes similar to the one achieved in a previous study (see Ref. Q). 

We conclude that accurate control of flow rates is crucial to achieving fast mixing. We rec¬ 
ommend that flow rates be analyzed by experimental quantitation of inlet velocities using, for 
example, micron-resolution particle image velocimetry (as performed by Hertzog et al. in Ref. Il2l ). 


C. Thermophysical Parameters 


We next studied the stability of the mixing time of the optimized mixer to changes in the 
thermophysical coefficients of the denaturant solution or in the concentration values needed to 
control the folding process. In physical experiments, these changes may result from uncertainties 


in conditions or solution properties (temperature, pressure, dilution, etc.) [2^ the following sections 


summarize. 


1. Denaturant Solution Parameters 


We chose for our work guanidine hydrochloride (GdCl) as a typical denaturant [l0l.l2ll] described 
by the following parameters: diffusivity in background buffer of D = 2 x 10“® m^ s“^, denaturant 
solution dynamic viscosity of 77 = 9.8 x 10“^ kg m“^ s“^ and mass density of p = 1010 kg m“^. 
The thermophysical properties of GdCl solutions vary with concentration and ambient temper¬ 


ature: consistent with the experimental work of Refs. 


24 I . I 25 I and 


2a, (i) the density of the GdGl 


solutions can vary within [1000,1700] kg m (ii) its viscosity can vary in [4,11] x 10 kg m ^ s 


1-1 0-1. 


and (iii) its diffusivity can vary in [1.9,13] x 10 ® m^ s ^ 


We considered the impact of these parameter variations on mixing time. We varied each par am- 
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TABLE III. Variation in percent of the mixing time value (relative to the value of the optimal mixer indicated 
by -) as a function of values of the side injection velocity Ug (m/s) and the injection ratio Uc/ug. 


ratio 

Uc (m/s) 

0.0097 0.0195 0.0292 0.0389 0.0973 

0.5 

26317 5444.8 1542.7 

668.7 

126.6 

1 

6136.1 

1518.6 

278.4 

134.2 

40.5 

1.5 

2776.7 

667.4 

109 

56.8 

20.8 

2 

1631.1 

388.7 

57.4 

30.5 

12.4 

2.5 

1083.7 

229.5 

34.6 

18 

7.8 

3 

775.2 

142.6 

22.2 

11.1 

4.9 

3.5 

583.3 

91.1 

14.8 

6.6 

3 

4 

456.8 

59.7 

9.8 

3.6 

1.5 

4.5 

368.9 

41.6 

6.2 

1.8 

0.4 

5 

304.6 

30 

3.6 

- 

0.5 

5.5 

256.6 

22 

1.6 

1.4 

1.3 

6 

219.5 

16.4 

1 

2.4 

1.9 

6.5 

190.4 

12.2 

1.1 

3.2 

2.4 

7 

166.8 

9 

2.1 

3.9 

2.9 

7.5 

147.5 

6.4 

2.9 

4.5 

3.3 

8 

131.6 

4.4 

3.6 

5 

3.6 

8.5 

118.1 

2.8 

4.2 

5.4 

3.9 

9 

106.7 

5 

4.7 

5.7 

4.2 

9.5 

97 

9.2 

9.4 

10 

42 


eter within the aforementioned intervals using seven equispaced values. All possible configurations 
of parameters values were studied, which represents a total of 343 evaluations of the model. Then, 
similar to the work presented in Section II V A 31 we computed the mean evolution of mixing time 
for each parameter value format both its lower and upper bound, and while varying the remaining 
coefficients to all their possible values. 

The variations of the mean mixing time of the diffusion, density and viscosity are presented in 
Figure [11] We see that the diffusion was the most sensitive parameter, and can increase mixing 
time by up to 0.3/rs. The other two coefficients maintained the mean mixing time close to O.l/rs. 
We note all of these values reasonable for the design as the order of mixing time is preserved. We 
further note increasing viscosity and decreasing diffusivity and density result in lower mixing time. 









22 



FIG. 10. Percent variation of the mixing time relative to the optimal mixer (represented by X) as a function 
of variations of the side injection velocity Us (m/s) and the injection ratio Uc/us- Precise control of flow 
rates are crucial in mixing experiments. 


The effect of decreasing diffusivity may at first seem counterintuitive, but mixing time is the result 
of a geometry- and flow-rate-dependent convective diffusion process. For example, high diffusivity 
can result in significant decreases of denaturant concentration within the early-focusing region of 
the center jet, where fluid velocities are still too low to stretch material interfaces and decrease 
diffusion lengths of the center jet. The latter effect is discussed by Hertzog et al. (2004) (e.g., see 
Figure 2 of that reference). 
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FIG. 11. Mean mixing time (in s) obtained as a function of denaturant solution a) diffusivity, b) mass 
density and c) dynamic viscosity, and fixing the other two parameters. 

2. Concentration threshold 

Finally, we characterize the sensitivity of the mixer to the maximum and minimum denaturant 
concentration values of our mixing time (see (4)). The original mixer was designed to trigger 
unfolding for a concentration reduction of 60%. We here consider mixing times for denaturant 
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concentration reductions ranging from 10% and 92%. To this end, for a particular threshold value 
denoted by 7 , we identify and Wy, such that = 'y and they produce the minimum mixing 

time value 


M^) 



d^/ 


( 8 ) 


where Ca^ and denote the Y-coordinates of the points situated along the streamline defined by 
the intersection of the two symmetry planes z = 0 fj,m and x = 0 ym, i.e. the y-axis, where the 
denaturant normalized concentration is and respectively. 

Results are presented in Figure [T 2 j The mixer exhibited mixing times lower than 0.4/is for up 
to a reduction of 90%. We conclude that it is a robust design as the maximum reduction allowed 
by the flow rate ratios in this mixer was 92%. For a 70% denaturant concentration reduction, we 
observed a O.l/rs mixing time. The latter can be compared to the mixer of Ref. [iT] which showed 
mixing times of l//s for the same denaturant concentration reduction [l^. 



FIG. 12. Mixing time (in s) as a function of the percentage reduction in concentration as defined in Section 
IIV C 21 The concentration values inherent to the definition of mixing have significant effect on mixing time. 
Note the maximum concentration reduction allowed by complete mixing far downstream is 92%. 
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V. CONCLUSIONS 


We presented a detailed study of the robustness and performance of a microfluidic mixer design 
first presented in Ref. Il7|. The mixer is for protein folding dynamics studies and can be used to 
initiate the folding process of a protein by diluting a local denaturant concentration in a short 
time interval. In Ref. [ITI, we used a 3D numerical model and showed the ideal mixer shape, flow 
control parameters, and expected thermophysical parameters, which resulted in a mixing time of 
0.10 /US. We here studied the robustness of this mixer relative to expected variations of these major 
design features. In particular, we studied (i) the uniformity of the mixing time through the center 
inlet flow and (ii) the sensitivity of the mixing time with respect to key mixer parameters. The 
uniformity study showed that mixing time is quite stable throughout the majority of the inlet 
stream (up to a distance from the walls of 0.4 ^m). With respect to design robustness, we found 
that the details of the mixer design in the region near the channel intersections are essential to the 
performance, i.e., the shape the minimum channel widths near this inlet, the inlet flow velocity 
ratio, and possible (unwanted) asymmetries in the fabrication. Other factors such as inlet channel 
angles, mixer depth (above a certain minimum), fluid properties, and denaturant concentration 
thresholds for protein folding have significantly weaker effect on mixing time. 

Our analyses may provide a guide to designers and fabricators of protein folding mixer devices, 
and can be used to evaluate trade-offs between manufacturing quality, precision of flow control, 
and expected performance. Our work also serves as a case study associated with the general design 
and performance prediction of microfluidic devices, and may serve as a guide to designing complex 
and optimal fluidic systems. In the least, the work highlights the complexity and importance of 
predicting and managing uncertainty in the performance of microfluidic systems. 

In Table HVl for each parameter, we provide the mean, maximum and standard deviation values 
of the mixing time percentage variation regarding the optimal mixer obtained with this sensitivity 
analysis. 
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TABLE IV. Summary of major findings of our sensitivity analysis: Mean (Mean), Standard Deviation 
(Dev) and worst-case Maximum (Max) values of the mixing time in percentage of base design. For the 
sake of completeness, we also report the optimal value of each parameter (Opt) as well as the range of the 
considered values (Range). 


Parameter 

Opt 

Range 

Mean 

Dev 

M2LX 

Intersection Angle 

e 

jh 

[0,2^/5] 

3 

3 

16 

Channel Intersection 


0.8 

[0,1] 

21 

17 

51 


0.7 

[0,1] 

30 

26 

79 

Channel Width 

Wc (Aim) 

2 

[1,4] 

50 

67 

150 

Ws (/iin) 

3 

[1,4] 

53 

57 

181 

Wo (Aim) 

10 

[2,18] 

28 

41 

101 

Mixer Depth 

Depth (Atm) 

10 

[8,12] 

13 

13 

34 

Symmetry 

Symmetry (%) 

0 

[0.5,50] 

216 

196 

542 

Injection velocities 


0.2 

[0.005,0.92] 

68 

133 

304 


5.2 

[0.5,9.5] 

51 

153 

669 

Physical coefficients 

D (m^ s ^) 

2 xlO"® 

[1.9,13] X 10-9 

155 

198 

307 

u (kg m^ s“^) 

9.8 xlO-4 

[4,11] X 10-4 

4 

4 

11 

p (kg m-3) 

1010 

[1000,1700] 

2 

2 

6 

7(%) 

60 

[10,92] 

129 

271 

1282 
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